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Abstract 


As the main contribution of this article, we establish an option on a credit 
spread under a stochastic interest rate. The intense volatilities in financial 
markets cause interest rates to change greatly; thus, we consider a jump 
term in addition to a diffusion term in our interest rate model. However, 
this decision leads us to a partial integral differential equation. Since the 
integral part might bring some difficulties, we put forward a fairly new nu- 
merical scheme based on the alternating direction implicit method. In the 
remainder of the article, we discuss consistency, stability, and convergence 
of the proposed approach. As the final step, with the help of the MATLAB 
program, we provide numerical results of implementing our method on the 
governing equation. 
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1 Introduction 


Spread options have great importance in financial markets. Modern and 
specialized investment in contemporary markets is based on this reality that 
we can represent strategies and models to invest in particular markets such 
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as the oil and exchange markets. These models are widely used in societies 
in which the inflation rate can change heavily during a short period, in short 
periods, and investors face less loss in these markets. 

For example, we can consider an investment in the oil market. If the 
investor defines their call (buy) and put (sell) positions according to the 
oil assets and their derivatives, they do not face many volatility challenges. 
Because the asset price in the oil market changes over a period, the price of 
another asset will change equally in the same direction with a little increase 
or decrement. Another example of such markets, occurring in some countries, 
is intense volatilities of exchange. In this case, if investors investment to buy 
and sell currency (e.g., uses the difference between Dollar and Euro), they 
will not face intense price volatilities of this market from other markets. 

The credit derivatives market has come under the spotlight recently. 
Credit derivatives pricing and valuation have received much attention over 
the last few years, and they are considered as one of the most significant new 
financial tools. The value of the portfolio with credit derivatives is extremely 
sensitive to changes in the spread between default-free bonds and default- 
able bonds. Hence the investors can minimize and manage the exposure to 
credit risk. Credit default swaps and credit spread options are the most 
popular form of credit derivatives. Credit spread options payment depends 
on particular credit spread or credit-sensitive asset prices. The default or 
credit rating downgrade risks can be compensated by corporate bond credit 
spread; it is an additional remuneration paid to investors. Credit spread 
can protect the holders of corporate bonds from the loss of the rising yield 
spread. Deng, Johnson, and Sogomonian [1] studied the spread options in 
the energy market. In 2015, Zhou et al. [2] suggested a credit spread option 
pricing model that depends on the time, interest rate, and the logarithm 
of the credit spread, which was solved by GARCH and LongstaffSchwartz 
methods. In 1995, Longstaff and Schwartz [3] studied that the credit spread 
logarithm follows the mean-reversion process and priced credit spread options 
on the assumption that the change in the distribution logarithm is normal. 
In 2012, Su and Wang [4] hypothesized that the interest rate would follow 
the Vasichek model and that the default intensity would follow the jump 
propagation trend. Tchuindjo [5], in 2011, proposed closed-form solutions for 
pricing credit risky discount bonds and their European call and put options 
in the intensity-based reduced-form framework, assuming the stochastic dy- 
namics of both the risk-free interest rate and the credit spread are driven by 
two correlated Ho-Lee models. 

These studies have not examined the abrupt changes in jumps and ex- 
changes. In this study, a model investigates based on abrupt changes in the 
exchange. Also, the ADI numerical method with expected convergence is 
applied. The innovation of this article is in proving the sustainability of the 
method. Finally, implementing the model on the real data is one of the most 
significant tasks in this article. The ADI method is an appropriate method 
for solving large matrix equations and numerically solving parabolic and el- 
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liptic differential equations. This method is applied to prevent solving a large 
linear system. The advantages of this method are the optimal accuracy of 
stability, low computational volume, and low computer memory requirement. 


This article is organized as follows: In section 2, the model, equation, and 
initial and boundary conditions will be introduced. In section 3, the ADI 
numerical method and its implementation on the model will be presented. 
Consistency, stability, and convergence of the method will be addressed in 
section 4. In Section 5, the proposed approach will be implemented on the 
real data. The last section provides some research titles for the followers. 


2 Modeling of the credit spread option pricing 


In this article, we consider the option pricing, according to the model pre- 
sented by Longstaff and Schwartz [3]. That is if we take X; as a credit 
spread option in a stock market, then its logarithm will satisfy the following 
stochastic differential equation [2]: 


dX = (a— bX)dt+ sdW,, 
also according to the Vasicek model, we have 
dr = (a — Br)dt + odW2, 


where the parameters a,b, s,a, 2, and o are constant numbers and W, and W2 
represent the Brownian motions and they are are independent variables with 
corr(W1, W2) = p. In mathematical modeling of financial topics, neglection 
of the jumps may lead to inaccurate results, by considering the economic and 
political problems, and also natural disasters such as earthquakes and storms 
as the main source of sudden movements, namely, jumps. So, in order to be 
able to present a pricing model for options, let us add the jump phrase to 
the Vasicek interest rate model. According to the information, we have 


dX = (a—bX)dt + sdWi, 


Mt 
dr = (a— Br)dt + odW2 + d(S— Zi), 
i=1 


where M; is a Poisson process with the positive density rate A that will be 
appeared soon, and the variables 7; represent the jumps and are independent 
and identically distributed. The credit spread option price demonstrated by 
P(X,r,t), equals to the solution of the following partial integral differential 
equation [6]: 
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1 1 
P, + (a — bX — qs)Px + (a — Br — qo)P, + posPxy + 5° Pax + 57 Por 
+00 


=(r+ AVP + P(X,r+y,t)f(y)dy = 0, 
0 


where f(y) is the Poisson distribution probability density function and y is 
the jump value. 

By applying the change of variable t = T — 7, for t € [0,7], the last equation 
changes into the following form: 


1 1 
P, =(a — bX — qs)Px +(a— Br—qo4 Pr posPxr + 58° Pxx 


ya) Per (r +A—1)P. (1) 


The initial and boundary conditions for 0 < X < K,0 < r < H, and 
0<-+t<T are as follows: 


0?P 0?P 

7x2 Oe" 7) a 0, 7x2 nT) = 0, 
0?P o?P 
pa (%9,7) = 0, pa (& HT) = 0, 
OP 


(K,7r,7) =0 OP (x Aw) =0, 


OX ” Or 


P(0,r,7) = P(X,0,7) = 1, P(X,7,0) = G(a), 


where G(x) = max{X — K,r — K,0} is the payoff function of an European 
put option (in Figure 1, we display payoff function based on data in Table 1) 
and Kk is the strike price. 


3 Numerical Solution 


In order to solve (1), we apply the ADI method. To do this, we define the 
operator LP as follows: 


LP=I*P+L"P4+1I*"P+L°*P+L"P+6. 


We define the last-mentioned components corresponding to the ones in (1) 
(see [7, 8]), 
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L* P = (a— bX — qs)Px, (2) 
LP = (a— Br — qo)P,, (3) 
L*" P = (p08)Pxy, (4) 
L*X p = 5(2)Pxx, (5) 
LP = (5(0)? + 5g) Pr (6) 
&=(r+A-1P, (7) 


while q is the market price of the risk. Approximation of these derivatives 
by finite differences results in [9, 16] as follows: 


1 1 1 
pts _ pn re p's 
(2): a 1 — (a — ba; — qs) +4 7 J 
T 
n+4 n+ nti 
=P; ° — Pi = ACP? = Ls *); 


(a — bX — qs)(Ar) 


where A = A . By rearranging the last equation, one has 
1 1 
(1+ A\P”* - APS = Pi, i =0,1,2,...,n, and the following system is 
obtained: 
qi 1 
(20% (1+ .A)Py;'* — API, ® = PR, 
e n+4 n+4 n 
i=l: (1+ A)P,; — AP, ° =Pi,, 
41 
Laos (1+ .A)P,;"* — AP; ® = PE, 
aed 1 
i=n—-1: (1+ A)Pn tj Ag Py 
i=n: (1+ A)Prt® — AP's = pr. 


The second term of the left in the last equation should be determined by 
known values as follows: 


+6 n+ n+ 
Poe -2P. t+ 1 1 i 
n+1jy nj n—-1j = n+ & =3 n+ se n+ ¢ 
is 05) 7-2 oUF 


n—l1j° 


Substitution the recent formula in the last equation (for 7 = n) results 


é=n:(1—A)POe 4 pers — p™ 


In the matrix representation, one gets 
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pet 
144 -A 0 0 0 a Poi 
0 1+A—-A 0 0 Pi; ; Pry 
0 0 14+A-A 0 a 7 PY 
0 0 0 1+A —A | | prtd PRs, 
0 0 0 A 1-A ne jee 
Pe 
Now we apply the finite difference for (3) as follows 
n+4 n+ 
Boot Pe 1 1 3 
ig ig r prt r prt 
— P. 6 P.. 3 
Ar a! Gg TSS ) 
n+ n+4 n+i n+ 
Ley ptt ag ° 1 3B atthe ij ) 
a k 7 k 
’ n+% ati n+i romti n+ rom 
= (1+3B)P,°° — 3B;Pi 3 = Pi’? — B,P§ =P, * — BPii8, 
where a = Fi*t |B; =a— fr; — qs, and fi; = prté — eee 
For j = 0,1, 2,...,m, we obtain the following system: 
pets "pnts 
n+ ‘4 n+4 
g=0: — (14+3B))Pig * — 38BoPiy* = foo, 
n 1 L & n+i 
gal: (1+3B))Py** -— 3B, PQ”! = fa, 
‘ n+ 1 pnt+4 n 
> j= es (1+ 3B,)P), * — 3BoPis* = fia, 
. ’ n+t ’ n+i 
j=m-1: (1+ 3B,,_ Pim =1 ~ 8Bin-1Pim * = fim-1, 
j=m: G4+an P08 9B fee 


By applying the boundary conditions, the second term of the left in the last 
equation (for 7 =m) should be determined by known values: 


al 
pee a opnts pes 


im+l1 im—l1°* 
Substitution recent formula in the last equation (for 7 = m) results in: 


Pp” 


ns ~ 3B, pnts = pets) = fim 


am 


j=m:(14+3B,,)P' 


=> (1-3B,, ieee a +3B,,Ps, im— = = fim- 
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This system can be written as the following: 


, , prts 
1+3B) -3By 0 0 se 0 ae Sio 
0 1+3B, -3B, 0 a 0 Fa * fa 
0 Oy TSR 6B, 3. 0 Prete fiz 
0 0 be. SG BASE ae 4 pe fim—1 
0 0 Qos 3B. 138, pret fim 
nj 


Using the finite difference method and the boundary conditions, (4) will 
transform into the following form: 


n+% 


n+4 
peta _ pe 
a ay Xr Xr p Xr pnts 
pipes 16 Lo re 
Ar = 75 * ) 
= (1-460) Pn*? — 230PMt, + baa +2304? +930"? 


+ ssp} — 230?" = fis, 


i i 2 H L d, posAt 
in which f;; are the sentences in n+ 3 and n+ g steps and C = 43,;°. Note 


that the finite difference in (4) follows the following formula: 


n+% n+ n+ n+s n+% n+4 
Pose Pp = Fey Eg Peg Sey Fy 
r hk . 
For 1 = 0,1,2,...,n, the following matrix equation is obtained: 
pe 
OF 
1 0 0 ee ee a 0 pete 
23C 1—46C 23C O ree eee tee vee 0 pete 
0 930 = 1—Ab6C 23C .-- +e ee dase 0 23 
DG S Fee, PM 
0 0 0 0 ++. 0 2301-460 23C] | Pn—2) 
0 0 0 0 --- 0 O 0 1 Pee 
we 
+3 3 
a 
230 23C ++ + oy | Po! 
~23C 230 «++ + 0 Pyj-i 
0 -—23C 23C --- 0 : 
i : 
= - BC 0 pee 
O D> ee ei 93G 936) | ented 
0 0 (0 «=930930) | Prt 


Mohamadinejad, Neisy and Biazar 


202 
n+4 
20 —23C” -0 0 0 pe fos 
dG 230.330 =e 0 tl As 
0 QO GC 930 2s. 0 ale 7 
+ 7 ; oo | a | 
0 0 0 23C —23C 0 uae a is 
0 0 0 O 20 —230 a fas 
0 0 0 0 23C —23C n—1j+1 ee 
prta 
nj+ 


Applying the finite difference for (5), we get the following system of equations, 
4=0,1,2,...,n 


n+2 n+s 
Pi 3 —P,; 2 
AT 
1 n n n+i n+2 
= 54 (-9E** Pi; +6 4 37LX% pnts _ 591 XX prt? 4 559E*% pts) 
Pree) 
cae Pog * = fog, i +2 +3 
Fea (1+ MOD) Pi"! — 55Po;"* — 5BP2;"* = fig, 
j=2 (1+110D)P> Pts —55P/,* — spits = fo;, 
=> 
pets n+2 nt+2 
jsm—2: (14 MOD) PAT y — SEPA TS — BBPR TT = fay. 
j=m-1: (1+ 10D)P, pred =55P 3 -55Py * = fn-iy; 
- n+2 
jon: PRYOR f,,, 


where fj; are all the sentences in the steps n+4 = ,n+3 = ,nt+é =,and D= ae 


The finite difference in (5) follows the following formula: 


Prag — 2Pig — Pity 
Pxx = el rom 25 


Write the finite difference, regarding (6), and then, we obtain the following 
results for 7 = 0,1,2,...,m: 
n+2 n+2 
Pi, °° — Pi? 1 n n n 
a (bP PB ec 127A, ame 2616L"" P,. +3 


Ar ~ 720 
—ar7anr pet + 1901n Pet), 


The resulted matrix, by substitution 7 = 0,1,2,...,m is as follows (EF = 


POOR (50° tr aE). 
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1 0 0 sie si : 0 
—1901F 1+ 3802E —1901FE 0 oo vee 0 
0 —1901E 1+ 3802E —1901E vee 0 
0 0 vee —1901E 1+ 3802E —1901E 0 
0 0 vee vee —1901F 1+ 3802E —1901E 
0 0 oe chee 0 0 1 
a 
prti fio 
ss fia 
Po 8 fia 
«for fal: 
pe fim-2 
eas f; - 
Pim—1 fi 
Bee 


jp (r+ -1)P 


5 
=> Prt = —(r +\—1ArPR + Pre. 


4 Consistency, stability, and convergence 


By proving the consistency and stability of the method, the convergence will 
be guaranteed via the Lax equivalence theorem [2]. 


4.1 Proof of consistency 


Let ® be the exact solution of (1), which depends on the independent vari- 
ables r, X, and t (e.g., L(®) = 0). Assume that the exact solution of the finite 
difference equation is ¢ (to write F(y) = 0) and that v is an arbitrary contin- 
uous function of r, X, and ¢t with a sufficient number of continuous derivatives 
such that L(v) can be evaluated at (in, je,kh). Then the error function at 
the point vi5,4 = v(in, je, hk) is defined as 75,.4(v) = F(vij,n) — L(vi,j,n). If 
Ti,j,n(U) + 0 as 7, € , and A tend to zero, then the finite difference equation 
is consistent with the partial-differential equation as follows: 
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Pi,j,kt+1 — i,j,k 


h 
i+1,j,k 7 Pi,j,k 1) Pij+1b 7 Pi,j,k 
_ (a bX; qs) Pi+1,5 Pig (a Br; go 4 va jt Pi,j 
n Xr € 
me Pit jt+Lk — Pitl1j,k — Pi,jtik + 201,54 — Pi-1,9,k — Pi,j—1,k + Pi-1,j-1,k 
Pp ne 
Le Pi+1,5,k — 20i,5,k + Pi-1,5,k 
2 "7 
1 1 Gig +Lk — 29i,5,k + Pij-1k 
(50? +5) pk TIEN (ry + A= Dpiiee (8) 


Let X; = Xo+in, let rj = ro+je, let A= a—bXo—qs, B= a+ }—Bro—4s, 
let C = pos, let D= 58°, let EF = 507 + xb, and let F = rp + A-—1. By 
considering the Taylor expansion and previous assumptions, (8) changes into 
the following equation: 


Op Wee . sn Pp | stop : 
aoe +t 6 ae + Ab gays + B-Beilgaa + 
2 93 2 93 
ir Oe a 
+ (A bnt) 7x3 I (B Bei) >a I ‘ 


By considering the smallest power of the above equation and assuming 
1 = € = kh (where « is a constant number), as h approaches zero, 7 will 
approach to zero too and the consistency of this method is proved. 


4.2 Proof of stability 


Proving stability directly from the definition is quite difficult, in general. 
Instead, it is easier to use tools from the Fourier analysis to evaluate the sta- 
bility of finite difference schemes. In this section, we analyze the stability of 
the ADI method. To study the stability, we use the Von Neumann approach 
[13, 14, 17]. The Von Neumann analysis is based on calculating the ampli- 
fication factor of schemes, G', and deriving conditions under which |G| < 1. 
We have introduced the numerical amplification factor G as 


m+1 
Ex 


G= ; 
ER 


(9) 


From (9), we may relate the error Ej’; at the nth time step with the initial 
error Ep ; [8] by 


n n 770 0 ika il 
Ex = GER, Bey = ere. 


For finding G, a simple procedure is to replace y;, in the scheme by 
G"e** ely for each k,l, and n, 
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n+1 _ prntl jika pily _ n 

VELY = Ger eet = GYR, 

n —_ an ik(a+An) ily __ pn jikAa ikea ily __ on (ikAa 
Pr+1 = Ge ( elt = Ge ee = py 1€ ; 

n —_ pin ike il(ytAy) _ pin_ika pily ilAy _.n WilA 
Yeiti = Greve (yt+Ay) — Greike pily pilAy _ euie (10) 


where Ax and Ay are positive. 

Without loss of generality of the problem, suppose that the coefficients of 
the relation (8) are assumed to be equivalent to constant A, where A is the 
largest value of the relation coefficients (8). 

Substituting (10) in relation (8) leads to the following relation: 


G-1 eikAx as | eilAy | eikAx eee ec tkAx eilAy a9 is etlAy 
— | \ \ 
At Ac  ' Ay | (Ax)? (Ay)? 
ei(kKAr+lAy) _ pikAx _ pilAy 12—e tkAx _ @—ilAy Le i(kAxv+lAy) 
t + (1 
tea 01) 


Suppose that kAx = [Ay = 9, that Ae = x = h, and that ae — any = 
n. Then by simplifying relation (11), G is obtained as follows: 


G=1+ A(—2h + 2hcos 6 — 2n(1 — cos20) + At) + 2Ahsin 01. 
We calculate G? and establish |G|? < 1 to obtain stability conditions from 
it (Because in this article, Az and Ay are considered values less than 1, we 
consider h < 7 in the |G|? formula) as follows: 
IG)? = (1 + A(—2n + 2n cos 6 — 2n(1 — cos20) + At))? + 4n? A? sin? 6. (12) 
If the following conditions are met, relation (12) will be stable [15, 5]: 


1+ A(—2n + 2n cos @ — 2n(1 — cos20) + At) < 2Ancos6, 4n?A? <1. (13) 


The stability conditions of the problem, according to relations (13), are as 
follows: 
At Zin 1 1 \ At Zuniat 1 1 
—, <min min 
Aa? — 2|A|’ A(2— Aw?)°? Ay? — 2|A|’ A(2 — Ay?) 


3 


Numerical results of this convergence and stability are given in the next 
section. 


5 Numerical results 


As an illustration, consider the credit spread and an interest rate to im- 
plement the ADI numerical method with 3 step-size. Then, we obtain the 
changes in credit spread and an interest rate for these markets. Table 1 shows 
the parameters used in this study. 
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Table 1: Data 


Parameters definitions values 
K Maxz(X) : maximum value of X 30 
H Max(r): maximum value of r 30 
qd the market price of the risk 1 
or volatility 0.3 
T Mazx(t): maximum value of t 0.5 
a,b, and s constant 0.5, 0.5, and 0.25 
r intensity rate of Poisson process 0.1 
p corr(W,, W2) 0.5 
r interest rate 0.03 


The payoff function can be illustrated in Figure 1 


Figure 1: Payoff function. 
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Pl 


[/Ne=30, Ny=20, e=0.90059 | 


$888 82 8 


f 


Be 


smo [Ne38, Nea. 0:3] 


NeXt, Npx30, 190 C087 


Figure 2: Numerical results using the ADI method with a step-size of é for n = 
0.3333, 0.5, and 0.66667 
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Zoom in 


Figure 3: Comparison of option prices for different interest rates such as r = 
0.03, 0.35, 0.55, and 0.9 


Eventually, regarding our model and the suggested numerical method, we 
can obtain the credit spread option price. As can be seen from the following 
figures, the price of the credit spread option varies for different amounts 
of the variable n (Figure 2 shows the credit spread option price functions 
at t = 0.33333, 0.5, 0.66667). Some option values for a different base on 
asset prices, for example, 70,110,170,220, and 270, are listed in Table 2 for 
different interest rates, for instance, r = 0.03, 0.35, 0.55, and 0.9. In addition, 
Figure 3 shows a comparison of option prices for the mentioned interest rates, 
which shows that interest rates and option prices are inversely related. 

By selecting Ax = Ay = 1 and At = i the error results of time steps were 
obtained 0.6515, 0.4798, 0.3173, 0.0213, and 0, respectively. As the time steps 
got smaller, these errors became smaller (by selecting Ax = Ay = 5 and At = 
a) the error results of time steps were obtained 0.3173, 0.0030, 0.0009, 0.0006, 
and 0, respectively), but in both cases, because our selected condition has 
been extracted from the stability range, the result of the error confirms the 
correctness of this range. As you can see, the resulting errors decrease and 
tend to zero. So, this method is convergent and stable [9]. 


Table 2: Option values for the different underlying asset prices 


P(x,r,7) x 

70 110 170 220 270 
r=0.03 925 1386.9 2079.8 2657.2 3234.7 
r=0.35 922.2 1382.8 2073.7 2649.4 3225.3 
r=0.55 920.6 1380.3 2069.9 26446 3219.3 
r=0.9 917.6 1375.9 2063.2 2636.1 3209 
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6 Conclusions 


The main goal of this article was to propose a pricing model, which is based 
on the credit spread and the interest rate, such that some jumps on the 
second are possible. Though the bringing jump in the modeling seems an 
interesting idea, it can bring some difficulties when it comes to implementing 
a numerical method. Therefore, to deal with this problem, as the second 
step, we postulate a fairly new numerical method call, ADI with step-size 
é combined with Adams—Bashforth formula [9]. To continue this work, one 
good idea is to use other methods such as the machine learning, Chebyshev 
cardinal wavelets, and wavelets Galerkin method to obtain the results [4, 
3]. Definitively, comparing the results of different methods could be highly 
instructive. 
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